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In  thermoelectric  generators,  the  heat  sources  are  usually  fluids  or  flames.  To  simplify  the  co-design  and 
co-optimization  of  the  fluid  or  combustion  system  and  the  thermoelectric  device,  which  are  crucial  for 
maximizing  the  system  performance,  a  three-dimensional  thermoelectric  generator  model  is  proposed 
and  implemented  in  a  computational  fluid  dynamics  (CFD)  simulation  environment  (FLUENT).  This  model 
of  the  thermoelectric  power  source  accounts  for  all  temperature  dependent  characteristics  of  the  mate¬ 
rials,  and  includes  nonlinear  fluid-thermal-electric  multi-physics  coupled  effects.  In  solid  regions,  the 
heat  conduction  equation  is  solved  with  ohmic  heating  and  thermoelectric  source  terms,  and  user  defined 
scalars  are  used  to  determine  the  electric  field  produced  by  the  Seebeck  potential  and  electric  current 
throughout  the  thermoelements.  The  current  is  solved  in  terms  of  the  load  value  using  user  defined  func¬ 
tions  but  not  a  prescribed  parameter,  and  thus  the  field-circuit  coupled  effect  is  included.  The  model  is 
validated  by  simulation  data  from  other  models  and  experimental  data  from  real  thermoelectric  devices. 
Within  the  common  CFD  simulator  FLUENT,  the  thermoelectric  model  can  be  connected  to  various  CFD 
models  of  heat  sources  as  a  continuum  domain  to  predict  and  optimize  the  system  performance. 

©  2010  Elsevier  Ltd.  All  rights  reserved. 


1.  Introduction 

Thermoelectric  generators  are  unique  power  sources  that  di¬ 
rectly  convert  heat  into  electricity  by  means  of  semiconductor 
materials.  They  are  of  great  interest  in  energy  applications  due  to 
their  well-known  merits  such  as  high  durability  and  environmen¬ 
tal  friendliness,  and  they  recover  thermal  energy  to  generate 
power  in  a  simple  manner.  Thermoelectric  devices  have  been  in¬ 
stalled  in  automobile  exhaust  pipes  to  reduce  the  power  load  on 
the  vehicle’s  alternator  [1],  and  in  biomass  or  gas  fired  heaters  to 
provide  power  for  fluorescent  lights,  TVs,  pumps,  fans  or  control 
panels  [2-6].  Another  kind  of  design  has  integrated  thermoelectric 
generators  in  various  micro-reactors  and  micro-combustors  as 
miniature  power  sources  for  portable/mobile  electronic  devices 
and  sensor  network  nodes  to  compete  with  low  energy  density 
electrochemical  batteries  and  fuel  cells  with  complex  on-board 
reformers  [7-15].  The  concept  of  combustion-driven  thermoelec¬ 
tric  generation  has  also  been  implemented  in  sophisticated  config¬ 
urations  of  active  counterflow  heat  exchangers  and  reciprocating 
flow  combustion  in  porous  materials  [16-18]. 

When  thermoelectric  devices  are  coupled  with  these  fluid  and 
combustion  systems,  the  interaction  between  the  heat  source 
and  the  generator  is  critical  to  the  overall  performance.  Although 
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the  power  generation  and  efficiency  of  thermoelectric  generators 
have  been  discussed  at  length,  their  influence  on  the  heat  source 
can  be  important  as  well  [1,3,9,12].  To  identify  the  impact  of  incor¬ 
porating  thermoelectric  generators  into  energy  systems  and  to 
design  a  new  generation  of  power  applications  with  increased  per¬ 
formance,  modeling  and  design  at  the  system  level  are  mandatory. 
In  practical  applications,  the  majority  of  heat  sources  for  thermo¬ 
electric  generators  are  in  fluid  form,  whereas  thermoelectric  effects 
are  described  by  solid  heat  transfer  terms.  Thus,  the  main  chal¬ 
lenge  in  the  development  of  such  system  models  is  to  simulate 
the  fluid-structure  coupled  multiphysics  effects,  and  the  tempera¬ 
ture  gradients  across  the  top  and  bottom  surfaces  of  the  thermo¬ 
electric  generator  and  across  the  contact  surfaces  of  the  hot  and 
cold  source  fluids  should  depend  upon  each  other,  i.e.,  using  the 
third  kind  of  thermal  boundary  conditions  instead  of  the  first  or 
second.  The  designer  must  integrate  the  thermoelectric  device 
model  into  the  fluid  model,  and  the  system  optimization  (combus¬ 
tion  control  and  combustor  design,  channel  and  heat  exchanger 
shape,  flow  rate,  inlet  direction,  etc.)  must  be  done  in  conjunction 
with  the  thermoelectric  generator  optimization  (geometry  of  pel¬ 
lets,  segmentation/cascade,  and  deployment  of  modules/thermo¬ 
couples,  among  others). 

Thus  far,  a  number  of  system  models  have  been  proposed  by 
coupling  the  analytical  thermoelectric  model  with  various  fluid 
models  obtained  from  simplified  assumptions  for  specific  applica¬ 
tions  [19-29].  These  system  models  are  useful  in  the  rough 
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Nomenclature 

C 

capacitance 

s 

uniform  cross-sectional  area  of  the  leg 

c 

specific  heat  of  the  element 

T 

temperature  distribution  of  the  leg 

I 

current 

t 

time 

j 

electric  current  density 

L 

temperature  of  the  cold  fluid  source 

k 

thermal  conductivity  of  the  material 

Tc 

temperature  of  the  cold  junction 

Ki 

interface  thermal  conductance  between  the  generator 

Th 

temperature  of  the  hot  junction 

and  the  hot  fluid  source 

Tw 

temperature  of  the  hot  fluid  source 

I<2 

interface  thermal  conductance  between  the  generator 

1/ 

electric  potential 

L 

and  the  cold  fluid  source 
length  of  the  p-  and  n-type  legs 

x,y,z 

directions  of  the  3D  coordinate  system 

m 

mass  density  of  the  element 

Greek  symbols 

P 

power  on  the  load 

a 

Seebeck  coefficient  of  the  device,  a  =  ap  -  an 

Q 

heat  flow  density  in  the  leg 

^p,n 

Seebeck  coefficient  of  the  p-  or  n-type  material 

Qc 

rate  of  heat  transfer  from  the  cold  junction  to  the  heat 

AT 

temperature  difference  across  the  element 

sink 

p 

efficiency 

Qh 

rate  of  heat  transfer  from  the  heat  source  to  the  hot 
junction 

P 

electrical  resistivity  of  the  material 

R 

generator  internal  resistance 

Subscripts 

Ri 

load  resistance 

P 

p-type 

Rp,n 

resistance  of  the  p-  or  n-type  leg 

n 

n-type 

estimation  of  interactions  between  the  fluid  sources  and  the  gen¬ 
erators  and  are  able  to  suggest  general  strategies  for  the  afore¬ 
mentioned  optimization.  However,  simplified  fluid  models  are 
not  sufficient  to  accurately  predict  the  detailed  behavior  of  most 
practical  fluid  systems  with  multidimensional  construction,  irreg¬ 
ular  geometry,  dynamic  variations  in  mass  flow,  or  complicated 
combustion  processes,  nor  can  they  precisely  transfer  nonuniform 
heat  flow  and  temperature  distributions  to  the  thermoelectric 
model  as  the  boundary  condition.  On  the  other  hand,  the  assump¬ 
tion  of  constant  material  properties  made  in  the  analytical  ther¬ 
moelectric  model  is  not  realistic  in  many  applications.  For 
thermal  fluid  tubes  and  fuel  fired  combustors  on  which  thermo¬ 
electric  generators  can  be  mounted,  the  temperature  therein  is 
high,  ranging  from  hundreds  to  more  than  one  thousand  degree. 
Thus  the  values  of  the  three  principal  properties  of  the  p- type 
and  n-type  materials,  i.e.,  aPi„,  pp>n,  and  kPtTl,  will  strongly  vary 
with  temperature. 

The  nonlinearity  in  thermoelectric  device  modeling  due  to  the 
temperature  dependency  of  material  properties  in  most  cases 
necessitates  a  numerical  approximation  instead  of  analytical 
methods,  whereas  complicated  fluid  power  systems  can  be  simu¬ 
lated  using  available  computational  fluid  dynamics  (CFD)  tech¬ 
niques.  Particularly,  a  number  of  CFD  models  have  been 
developed  [15,30-34]  to  deal  with  the  thermoelectric  generation 
by  FLUENT,  a  widely  used  industry  code.  However,  these  CFD  mod¬ 
els  only  simulate  the  status  of  heat  source  fluids  or  convection / 
radiation  effects  on  thermoelements  without  an  appropriate  model 
of  thermoelectric  generators.  In  other  recent  works,  both  FLUENT 
and  thermoelectric  generator  simulation  have  been  employed  to 
study  the  system  thermal  behavior,  but  coupling  of  the  two 
numerical  programs  has  not  been  completed,  i.e.,  the  interaction 
is  in  a  single  direction  [35,36].  Although  in  principle,  numerical 
heat  transfer  schemes  of  thermoelectricity,  such  as  the  model  in 
the  commercially  available  finite  element  method  (FEM)  program 
ANSYS  [37]  (used  in  [35]),  the  one  used  in  [36],  and  those  described 
in  [38-40],  can  be  linked  to  a  CFD  simulator  through  a  relaxation 
strategy,  auxiliary  coding  is  usually  required  for  iteration  control, 
data  transfer,  and  synchronization  between  different  simulators. 
The  processing  of  coupling  relationships  of  multidimensional 
boundary  values  for  the  same  model  domain  but  with  different 
simulators  may  be  especially  difficult  [41  ].  Therefore,  it  is  certainly 


preferable  to  model  both  fluid  behaviors  and  thermoelectric  power 
output  with  the  same  tool  to  minimize  the  number  of  simulators 
with  high  usability  and  without  loss  of  accuracy. 

This  work  presents  a  three-dimensional  (3D)  numerical  solu¬ 
tion  to  the  fluid-structure  coupled  problem  by  implementing  a 
FLUENT  compatible  thermoelectric  model.  FLUENT  and  its  user-de¬ 
fined  functions  (UDFs)  were  selected  for  fuel  cells  to  model  the 
electrochemical  reactions  and  the  electric  field  caused  by  the 
Nernst  potential  throughout  the  cell  [42].  The  most  significant 
benefit  of  using  FLUENT  to  model  such  power  sources  is  that  the 
existing  knowledge/experience  of  this  computational  tool  and 
the  innumerous  fluid  models  already  implemented  by  it  can  be 
efficiently  utilized  to  predict  the  coupled  system  performance. 
Moreover,  UDFs  allow  for  a  convenient  customization  and 
enhancement  of  FLUENT,  where  the  thermoelectric  device  model 
can  be  connected  to  the  thermal  fluid  model  as  a  continuum  do¬ 
main  to  avoid  boundary-value  transport  problems,  making  the 
co-design  of  the  thermoelectric  generator  and  fluid  heat  sources 
simpler  and  less  time-consuming.  In  the  following  sections  of  this 
paper,  3D  governing  equations  and  the  general  multidimensional 
numerical  algorithm  of  thermoelectric  generation,  model  imple¬ 
mentation  in  FLUENT,  and  model  results,  will  be  described. 


2.  Multidimensional  multiphysics  numerical  scheme 

A  typical  thermoelectric  power  system  is  shown  in  Fig.  A.l,  in 
which  many  n-  and  p- type  semiconductor  legs  composing  the 
generator  are  connected  thermally  in  parallel  between  the  hot 
and  cold  fluid  sources  and  electrically  in  series  to  power  the  load 
circuit.  Details  of  the  general  numerical  algorithm  used  to  deal 
with  nonlinear  issues  such  as  heat  production  due  to  the  Joule  ef¬ 
fect,  unsteady  state  conduction,  and  temperature  dependent 
material  properties  in  one-dimensional  (ID)  transport  equations 
for  such  thermoelectric  generators  were  documented  in  a  previ¬ 
ous  study  [43],  in  which  we  introduced  a  combined  finite  differ¬ 
ence  and  Newton-Raphson  method  based  on  the  governing 
equation  system  from  the  energy  conservation  theory.  The  main 
intention  of  Section  2  is  to  extend  the  ID  algorithm  to  the  multi¬ 
dimensional  implementation  of  the  multiphysics  numerical 
scheme  for  thermoelectric  generators  operating  in  thermal  fluids. 
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As  is  shown  in  Fig.  A.l,  it  is  clear  that,  for  fluid  heat  sources  with 
flowing  mass,  neither  their  temperature  nor  their  heat  flow  profiles 
along  the  contact  interface  to  thermoelectric  generators  can  be 
uniform.  Obviously  the  results  of  the  thermal  analysis  are  now 
multidimensional,  e.g.,  the  temperature  and  heat  flow  distribution 
should  be  expressed  as  Tpn(x,y,z)  and  gpn(x,y,z)  in  a  3D  model, 
respectively.  In  the  case  of  temperature  dependent  material  prop¬ 
erties,  neglecting  the  side  heat  loss  outside  the  solid  legs,  the  3D 
governing  equation  of  the  thermal  field  inside  a  control  volume 
is  written  as, 

dT 

^Qp,n  =  —  V[kp,n(Tp,n)VTpji]  =  TYlpnCpn  ^  +  J  Ppn(JP,n) 

+  [Va  p,n(Tp,n)]W,  (1) 

where  the  first  term  in  the  right  side  is  the  transient  term,  the  sec¬ 
ond  term  is  the  temperature  dependent  Joule  electrical  energy  pro¬ 
duction,  and  the  third  source  term  represents  the  contributions 
from  Peltier  (inhomogeneous  or  segmented  materials)  and 


Thomson  effects.  By  an  electro-thermal  analogy,  the  heat  transfer 
governed  by  (1)  is  illustrated  in  Fig.  A.2  (a),  a  control  volume  of 
the  3D  thermal  resistor  network. 

The  electrical  analysis  requires  these  thermal  analysis  results, 
mainly  the  temperature  profiles,  to  find  the  total  Seebeck  voltage 
generated  and  the  temperature  dependent  material  properties. 
Keeping  the  non-ohmic  current-voltage  relation  in  both  legs  of 
the  device  in  mind  [44],  the  governing  equation  of  the  electric  field 
inside  a  control  volume  under  steady  state  is  written  as, 

VV  =  —0CPjn(TPjn)VTPjn  —  Ppji(TP:n)Ji  (2) 

where  the  first  term  in  the  right  side  of  (2)  is  the  Seebeck  electro¬ 
motive  force  (EMF)  increase  due  to  the  temperature  gradient,  and 
the  second  term  is  the  voltage  drop  due  to  the  current  flowing 
through  the  control  volume.  As  a  result  of  the  3D  temperature  dis¬ 
tribution,  the  electric  current  and  potential  distribution  are  also  3D, 
i.e.,7(x,y,z)  and  V(x,y,z).  In  addition,  the  values  of  Rp<n  must  be  ob¬ 
tained  by  a  calculation  of  a  3D  resistor  network  for  the  3D  resistiv¬ 
ities  pPin(x,y,z),  which  is  shown  in  Fig.  A.2  (b),  where  the  current 


Fig.  A.2.  A  control  volume  of  (a)  3D  thermal  resistor  network  Eq.  (1),  (b)  3D  electric  resistor  network  Eq.  (2). 
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vector  flows  in  three  directions  against  the  Seebeck  voltage  sources. 
Although  the  ID  analysis  is  fast  and  easy  to  manipulate,  it  lacks  de¬ 
tails  about  the  lateral  distributions  of  electric  potential  and  current, 
and  hence,  it  is  impossible  to  obtain  the  effective  electrical  resis¬ 
tance.  Therefore,  the  3D  electric  field  analysis  of  thermoelectric 
generators,  by  which  the  main  results  of  power  and  efficiency  can 
be  obtained,  should  be  accompanied  by  the  3D  thermal  field 
analysis. 

The  complete  solution  of  a  multidimensional  multiphysics 
problem  relies  on  a  clear  data  and  algorithm  flow  chart,  in  which 
various  numerical  approximation  methods  for  the  nonlinear  differ¬ 
ential  equation  system  can  be  identified  and  connected  as  a 
sequential  process.  Fig.  A.3  shows  such  a  self-consistent  flow  chart 
for  thermoelectric  generators  operating  in  the  system  mode  as 
shown  in  Fig.  A.l.  The  basic  mechanism  consists  of  four  simulation 
modules:  a  3D  CFD  simulation  is  used  for  the  energy  and  mass 
transfer  status  in  the  hot  and  cold  fluid  domains;  a  3D  thermal  sim¬ 
ulation  for  the  temperature  and  heat  conduction  distribution  in  the 
thermoelectric  domain;  a  3D  electric  simulation  for  the  electric 
field  distribution  and  electric  conduction  in  the  thermoelectric  do¬ 
main;  and  a  ID  circuit  simulation  for  the  current  and  power  perfor¬ 
mance  of  the  actual  load.  The  four  simulation  modules  are  coupled 
together  and  there  are  three  direct  coupling  relationships  among 
them.  First,  the  CFD  simulation  is  coupled  with  the  thermoelectric 
thermal  simulation  through  the  contact  boundary,  that  is,  the  3D 
boundary  temperature  and  heat  flux  of  both  the  fluid  domains 
and  the  solid  domain  should  be  kept  continuous  with  each  other 
as  a  whole.  Second,  within  the  thermoelectric  device,  the  thermal 
simulation  and  the  electric  simulation  are  coupled  with  each  other 
over  the  entire  3D  solid  domain  through  the  update  of  all  of  the 
temperature  dependent  material  properties,  the  temperature  dis¬ 
tribution,  and  the  current  density.  Third,  the  3D  electric  simulation 
is  coupled  with  the  load  circuit  as  a  field-circuit  coupling,  where 
the  overall  Seebeck  voltage  and  effective  internal  resistance,  calcu¬ 
lated  from  the  3D  electric  simulation,  are  applied  as  a  DC  voltage 
source  to  the  external  load  circuit. 

Unlike  the  thermoelectric  power  source,  the  load’s  current- 
voltage  relation  fulfills  Ohm’s  law,  and  there  is  no  need  to  consider 
the  power  output  in  a  3D  way.  Depending  on  the  position  of  the 
electrode  pads  of  the  actual  device,  the  translation  between  the 
3D  current  distribution  from  the  thermoelectric  electric  simulation 
and  the  ID  current  from  the  circuit  simulation  can  be  easily  carried 
out  because  the  current  is  the  only  boundary  condition  in  this 
field-circuit  coupling.  For  the  coupling  between  the  fluid  field 
and  the  thermoelectric  field,  however,  the  boundary  condition 
translation  is  much  more  difficult  because  the  CFD  simulation 
and  the  thermoelectric  thermal  simulation  both  output  3D  heat 
flow  and  temperature  distributions  at  the  boundary  surface.  There 


are  two  boundary  conditions  for  the  same  boundary  domain,  and 
one  of  them  must  be  known  beforehand  for  both  simulators  to 
start  the  thermal  simulations.  Fig.  A.3  shows  one  iteration  strategy, 
where  initial  values  of  the  boundary  heat  flow  distribution 
Qh,cix^y^z)  must  be  given  before  the  CFD  simulation  starts  to  run. 
Then,  the  temperature  distribution  Twa(x,y,z),  as  the  solution  of 
the  energy  transfer  equations  solved  in  the  CFD  simulation,  can 
be  used  in  the  numerical  thermal  analysis  of  the  thermoelectric 
field,  whose  analysis  results  are  fed  back  to  the  CFD  simulation 
to  update  Q/i,c(x,y>z)  as  the  new  boundary  values. 

In  spite  of  the  iteration  instance,  the  translation  of  the  multidi¬ 
mensional  thermal  boundary  condition  in  practical  simulation 
usually  involves  complicated  mathematics  and  computational 
technique,  as  shown  in  [41  ]  and  other  open  reports.  To  avoid  any 
theoretical  difficulties  brought  about  by  such  a  manipulation  be¬ 
tween  CFD  and  heat  transfer  simulators,  it  is  advantageous  to 
implement  the  thermoelectric  model  in  the  CFD  framework  as 
well;  thus,  the  solid  and  the  fluids  are  defined  in  a  continuum  zone 
and  the  boundary  is  simply  solved  as  interior  nodes  by  a  single  sol¬ 
ver.  The  next  section  will  present  the  detailed  treatment  of  ther¬ 
moelectric  modeling  in  FLUENT  to  enable  an  efficient  fluid- 
thermal-electric  coupled  simulation  as  the  scheme  displayed  in 
Fig.  A.3.  The  resistance  of  the  load  circuit  is  extracted  and  used  in 
the  thermoelectric  modeling  such  that  the  field-circuit  coupling 
is  also  included  to  some  extent. 


3.  Numerical  model 

The  3D  thermoelectric  generator  model  is  implemented  in  FLU¬ 
ENT  6.3,  a  finite  volume  method  (FVM)  CFD  package,  in  which  the 
solid  thermoelectric  phenomena  and  the  algorithm  described  in 
Fig.  A.3  are  taken  into  account  by  developing  UDF  through  ANSI- 
C  language.  The  model  serves  as  an  add  on  module  in  parallel  with 
other  power  source  modules  [42]  provided  within  the  standard 
FLUENT  software. 

Meshing  of  the  GAMBIT  pre-processor  of  FLUENT  is  applied  to 
leg  volumes  in  which  tetrahedral  elements  are  dominant  in  this 
work,  but  one  may  of  course  use  other  meshers  and  element  types 
if  the  device  under  study  has  irregular  leg  shapes.  The  grid  gener¬ 
ated  is  read  into  the  FLUENT  solver  and  scaled  to  align  the  geomet¬ 
rical  units.  Both  the  thermal  analysis  and  the  electric  analysis  use 
the  same  grid,  and  thus,  they  have  the  same  3D  coordinate  system 
with  regard  to  the  cells  and  boundary  walls  of  the  leg  volumes. 

FLUENT  has  the  ability  to  compute  the  conduction  of  heat 
through  solids  when  the  Energy  Model  is  activated.  Temperature 
dependent  thermal  conductivities  can  be  specified  by  polynomial 
functions  and  assigned  to  p-  and  n-leg  cell  zones  for  p-  and  n-type 
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Fig.  A.3.  Energy  conservation  based  multidimensional  multiphysics  simulation  procedure  for  thermoelectric  fluid  power  systems. 
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materials,  respectively.  The  thermal  boundary  condition  of  the 
thermoelectric  generator  does  not  need  to  be  specified  in  the  sys¬ 
tem  modeling  because  it  will  be  solved  in  terms  of  fluid  system 
boundary  conditions,  and  all  aspects  of  fluid  flow,  heat  and  mass 
transfer  in  heat  sources  are  handled  by  FLUENT.  The  transient 
source  term  in  (1)  can  also  be  included  automatically  with  appro¬ 
priate  settings  in  FLUENT.  The  values  of  Joule,  Thomson,  and  Pel- 
tier,  i.e.,  steady  state  source  terms,  however,  rely  on  the  3D 
electric  analysis  results,  as  shown  in  (1)  and  Fig.  A.3.  Thus  an  elec¬ 
tric  field  solution  is  required  to  complete  the  thermal  field 
calculation. 

FLUENT  solves  generic  transport  equations  for  user  defined  sca¬ 
lars  (UDS)  which  can  be  used  to  determine  the  electric  field  in  a 
CFD  environment.  The  difficulty  herein  is  that  the  governing  Eq. 
(2)  is  not  in  the  generic  transport  equation  form  due  to  the  Seebeck 
EMF.  To  correctly  write  (2)  for  the  solver,  two  UDS  are  used  to  rep¬ 
resent  the  two  terms  in  the  right  side  of  (2).  UDSO  is  used  to  rep¬ 
resent  the  3D  potential  distribution  of  the  ohmic  voltage  drop 
caused  by  the  current  vector, 

VUDSO  =  -Pp,n(Tp,n)J.  (3) 


conditions  for  UDS1  is  the  same  as  that  for  UDSO  except  that  a  con¬ 
stant  value  is  used  to  calculate  the  flux  J  instead  of  J.  Thus,  the  elec¬ 
tric  field  calculation  function  in  Fig.  A.3  can  be  implemented  with 
regard  to  UDS1  if  the  source  term  is  appropriately  included. 

Before  proceeding  to  the  implementation  of  the  right  side  of  (7), 
we  recall  the  source  term  in  the  3D  thermal  Eq.  (1).  An  energy_- 
source  UDF  is  written  to  modify  the  heat  conduction  computation 
to  include  the  Joule,  Thomson,  and  Peltier  heating.  FLUENT  calls  the 
UDF  as  it  performs  a  global  loop  on  cells  to  compute  the  source 
term  and  returns  it  to  the  solver.  Because  the  current  distribution 
calculation  has  been  done  with  UDSO,  the  Joule  source  term  can  be 
easily  obtained.  In  the  practical  implementation,  the  Joule  heat  can 
be  calculated  as  the  product  of  the  sum  of  the  squares  of  vector 
components  of  the  gradient  of  UDSO  and  the  temperature  depen¬ 
dent  electrical  conductivity, 


J  Pp,n(Tp,n) 


1 

Pp,n(Tp,n) 


(VUDSO)2 


1 

Pp,n  (Tp,n) 


fd  UDS0\2  ( <9UDS0\ 2 

\~W~) 


/  <9UDS0\ 2 

\^r) 


(8) 


Solving  for  3D  electrical  conduction  is  directly  analogous  to  the 
computation  of  heat  transfer.  The  field  throughout  the  conductive 
regions  is  calculated  based  on  flux  (charge)  conservation  in  each 
cell, 

V/  =  0,  (4) 

so  we  have 


The  source  terms  of  Thomson  and  Peltier  (in  the  case  of  inhomoge¬ 
neous  or  segmented  materials)  in  each  cell  involve  the  gradient  of 
the  Seebeck  coefficient  ap>„.  To  express  VaPi„(TPin)  for  the  source 
term,  UDS2  is  used  to  represent  the  scalar  of  aPi„(rPin).  We  consti¬ 
tute  a  transport  equation  for  UDS2, 

V(-VUDS2)  =  0,  (9) 


V 


1 

Pp,n(Tp,n) 


VUDSO 


=  0. 


(5) 


FLUENT  solves  this  Laplace  equation  for  the  potential  field  by 
enforcing  a  specific  flux  J  on  one  boundary  wall  of  both  p-  and  n- 
type  solid  regions  and  zero  potential  on  the  other  boundary  wall 
to  represent  the  ground,  where  the  I  value  of  the  present  iteration 
is  the  basis  on  which  the  specific  J  distribution  on  the  UDSO  bound¬ 
ary  is  calculated.  The  temperature  dependent  electrical  conductivi¬ 
ties  l/pp,n(TPin)  are  specified  as  the  diffusivities  of  UDSO  by 
polynomial  functions  for  p-  and  n-type  materials.  Referring  to 
Fig.  A.3,  the  functions  of  the  resistor  network  calculation  and  cur¬ 
rent  distribution  calculation  are  both  implemented  in  terms  of 
UDSO. 


where  a  unit  diffusivity  for  UDS2  is  taken.  It  should  be  remarked 
that  the  result  for  UDS2  from  the  solver  is  not  interesting  because 
(9)  does  not  hold  any  physical  meaning.  Instead,  once  per  iteration, 
an  at_end  UDF  is  executed  to  fill  UDS2  with  the  local  Seebeck  coef¬ 
ficient  aPin(Tp>n)  according  to  the  temperature  TPi„(x,y,z)  for  each 
cell.  The  solved  scalar  field  for  UDS2  is  replaced  therein,  where  tem¬ 
perature  dependent  Seebeck  coefficients  are  specified  in  the  UDF  by 
polynomial  functions  forp-  and  n-type  materials,  respectively.  If  we 
substitute  VUDS2  into  the  source  term  of  Thomson  and  Peltier  in 
(1),  we  obtain  the  thermoelectric  source  term, 

[Voe p,„(Tp,n)]TP,J  =  -VUDS2  VUDSO.  (10) 

Pp,n  \*  P,n) 

It  consists  of  three  components, 


UDS1  is  used  to  represent  the  3D  Seebeck  EMF  distribution  pro- 

Pttp,n(Tpji) 

II 

<9UDS2 

Tp,n 

<9UDS0 

(Ha) 

duced  by  the  temperature  field, 

dx 

dx 

Pp,n(Tp,n) 

dx  ’ 

-VUDS1  =  a :P,„(TP:„)VTPA, 

(6) 

9ttp,n(Tp,n) 

II 

<9UDS2 

TPin 

(9UDS0 

(lib) 

dy 

dy 

Pp.nTp.n) 

dy  ’ 

whose  divergence  is 

Pttp,n(Tpn) 

II 

<9UDS2 

TP,n 

<9UDS0 

(11c) 

V(-VUDSl)  =  V[ap,n(TPin)Vrp,n]. 

(7) 

dZ 

dz 

Pp,n(Tp,n) 

dz  ’ 

If  we  take  an  unit  diffusivity  for  UDS1,  (7)  becomes  a  generic  trans¬ 
port  equation  where  the  right  side  can  be  assumed  to  be  the  source 
term,  although  it  does  not  have  a  real  physical  meaning  as  a  source 
in  the  thermoelectric  modeling.  The  setting  of  the  boundary 


all  of  which  are  included  in  the  sum  of  the  energy_source  UDF. 

Now  we  proceed  to  the  source  term  in  (7)  for  UDS1,  imple¬ 
mented  by  an  UDS1  .source  UDF, 


Table  A.1 

Transport  equation  parameters  of  scalar  fields. 


Diffusivity 

Source  term 

Top  wall  boundary 

Bottom  wall  boundary 

UDSO. 

i 

Pp,„(Tp ,n) 

0 

specified  flux  (/) 

specified  potential  (0) 

UDS1. 

1 

UDSl_source  UDF 

specified  flux  (0) 

specified  potential  (0) 

UDS2. 

1 

0 

specified  flux  (0) 

specified  potential  (0) 

UDS3. 

1 

0 

specified  flux  (0) 

specified  potential  (0) 

UDS4. 

1 

0 

specified  flux  (0) 

specified  potential  (0) 

UDS5. 

1 

0 

specified  flux  (0) 

specified  potential  (0) 

Temperature. 

^p,n(Lp,n) 

energy_source  UDF 

from  CFD 

from  CFD 
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V[0Cp>n  (Tpn)  VTpn] 


dttp  n^Tpn)  dMp,n(Tp,n)  ^ 

dx  +  dy 

dOtp^n(Tp:n)  ~~§£~ 

dz 


respectively,  as  three  scalar  fields.  Similarly  to  UDS2,  we  constitute 
transport  equations  for  these  scalars, 


V(-VUDS3)  =  0, 
V(-VUDS4)  =  0, 
V(-VUDS5)  =  0, 


(13a) 

(13b) 

(13c) 


To  express  the  above  equation  in  the  UDF,  UDS3,  UDS4,  and  UDS5  not  for  their  solutions  from  FLUENT  but  to  fill  them  with  products  of 

are  used  to  represent  ocp,n(TPtn )^g5f  (xP,n(Tp,n) and  oc p,n(Tp,n)^-,  the  temperature  dependent  Seebeck  coefficient  and  the  three 


Table  A.2 

Comparison  of  ID  simulation  results,  Th  =  423  K,  Tc  =  303  K. 


Quantity 

Analytical 

ANSYS 

FLUENT  coarse  grid 

FLUENT  medium  grid 

FLUENT  refined  grid 

Measurement 

Qh,W 

81.3 

83.29 

80 

81.6 

81.8 

70 

P,  W 

3.98 

3.76 

3.08 

3.57 

3.62 

2.51 

Yj% 

4.89 

4.51 

3.85 

4.38 

4.43 

3.6 

I,  A 

1.08 

1.05 

0.952 

1.024 

1.032 

0.86 

Fig.  A.4.  Simulation  of  surface  heat  flux  (W/m2)  distributions.  Th  and  Tc  are  423  K  and  303  K,  respectively,  (a)  using  the  coarse  grid,  (b)  using  the  refined  grid. 
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components  of  the  temperature  gradient  in  the  at_end  UDF.  Be¬ 
cause  the  gradients  of  UDS3,  UDS4,  and  UDS5  can  be  accessed,  the 
source  term  (12)  for  UDS1  is  appropriately  implemented. 

At  the  end  of  every  iteration  (at_end  UDF),  both  the  ohmic  volt¬ 
age  drop  represented  by  UDSO  and  the  Seebeck  EMF  represented 
by  UDS1  are  determined.  Their  spatial  addition  in  terms  of  (2)  in 
the  full  solid  region  gives  the  3D  potential  distribution  V(x,y,z), 
and  the  effective  internal  resistance  of  the  legs,  Rp  and  Rn ,  can  also 
be  calculated.  Including  the  load  Ri ,  the  new  value  of  /  is 


_  UDSlp  -  UDS1 

Rp  +  Rn  +  R 


which  will  be  used  in  the  next  iteration  as  the  boundary  condition 
for  UDSO.  UDSlp  n  is  the  total  built-in  Seebeck  EMF  of  the  leg,  which 
should  be  equivalent  to  the  open  circuit  voltage  at  no  load.  This  cal¬ 
culation  is  relevant  to  the  position  of  the  actual  electrodes  at  which 
I  is  applied,  but  in  this  work,  the  UDS1  potential  of  the  boundary 
wall  is  simply  used.  The  previous  ID  study  [43]  aims  for  the  co-de¬ 
sign  in  which  the  thermoelectric  generator  system  is  associated 
with  complicated  electrical  systems.  For  the  3D  FLUENT  model,  if 
Ri  can  emulate  the  input  impedance  of  the  load  circuit,  the  afore¬ 
mentioned  field-circuit  coupling  is  also  included. 

As  far  as  can  be  observed,  all  thermoelectric  function  modules 
in  Fig.  A.3  are  implemented  by  UDF  and  UDS.  The  boundary  condi¬ 
tion,  diffusivity,  and  source  term  of  the  UDS  and  temperature  fields 
are  summarized  in  Table  A.l .  The  model  can  be  easily  implemented 
for  practical  devices  with  multiple  thermocouples.  When  the 
boundary  walls  of  the  legs  are  coupled  with  the  fluid  flows,  ther¬ 
mal  boundary  conditions  are  no  longer  required  on  them  because 
FLUENT  will  calculate  the  heat  transfer  directly  from  the  solution 


in  the  adjacent  cells  of  both  sides.  The  iteration  is  automatically 
done  by  FLUENT  until  convergence  is  achieved. 


4.  Model  results 

First,  a  comparison  of  the  ID  steady-state  simulation  between 
the  FVM  FLUENT  and  FEM  ANSYS  [37]  is  carried  out  to  theoretically 
validate  the  proposed  model.  The  example  considered  is  the  per¬ 
formance  of  the  thermoelectric  generator  described  in  the  previous 
study  [43],  the  TEC1-12706  thermoelectric  module  made  by  Tian¬ 
jin  Institute  of  Power  Sources,  China,  which  was  chosen  for  the 
experimental  validation  therein.  The  element  length  is  L  =  1.6 
mm  for  both  the  p- type  and  n-type  semiconductors,  and  the  ele¬ 
ment  cross-sectional  areas  are  Sp  =  Sn  =  1.4  x  1.4  mm2.  The  genera¬ 
tor  has  127  thermocouples  and  is  connected  to  a  linear  load 
resistance  Ri  =  3.4Q.  Second,  third,  or  fourth  order  polynomial 
functions  are  used  to  fit  the  temperature  dependency  of  the  p- type 
and  n-type  material  properties. 

Initially,  the  constant  temperature  boundary  condition  is  set  on 
the  top  and  bottom  surfaces  for  the  ID  simulation,  where  zero  heat 
loss  from  the  other  four  sides  of  the  legs  is  assumed.  With  the  tem¬ 
perature  polynomial  functions  used  for  varying  material  proper¬ 
ties,  the  FLUENT  thermoelectric  model  is  performed  to  analyze 
Qh ,  P,  ri ,  and  /  for  the  generator  operating  between  Tc  =  303  K  and 
Th  =  423  K.  A  gird  sensitivity  analysis  is  performed  for  the  pre¬ 
sented  FLUENT  model  to  check  its  gird  independency,  where  three 
grid  schemes  used  in  the  computation  are  designated  as  a  coarse 
grid,  a  medium  grid,  and  a  refined  grid,  respectively.  The  parame¬ 
ters  computed  by  FLUENT  with  the  three  grids  and  ANSYS  as  well 
as  the  analysis  using  the  material  properties  evaluated  at  an 


hot  junction/device  top  temperature  (K) 


Fig.  A.5.  Comparison  between  simulation  and  measured  results  for  various  temperature  difference.  Cold  junction/device  bottom  temperature  is  fixed  at  303  K,  Rt  =  3.4 Q.  (a) 
Qh,  (b)  P,  (c)  rj,  and  (d)  I. 
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average  temperature  of  363  K  are  summarized  in  Table  A.2.  The 
experimental  data  collected  in  [43]  are  also  provided  for  compari¬ 
son.  Clearly,  the  analytical  model  is  essentially  less  realistic  than 
the  numerical  models  and  yields  a  higher  r\  if  average  material 
properties  are  used  rather  than  the  temperature  dependent  mate¬ 
rial  properties.  This  result  has  been  pointed  out  in  [37],  where  the 
cause  of  the  efficiency  decrease  is  mainly  the  heat  evolution  of  the 
Thomson  effect.  When  the  grid  is  sufficiently  dense,  FLUENT  can 
output  almost  exactly  the  same  results  as  those  achieved  by 
ANSYS.  The  discrepancy  caused  by  the  relatively  coarse  grid  attri¬ 
butes  to  the  inherent  multidimensional  features  of  FLUENT.  In 
Fig.  A.4,  it  can  be  seen  that  the  nonuniformity  of  the  surface  heat 
transfer  contributing  to  the  interior  conduction  in  the  legs  of  a  sin¬ 
gle  couple  is  appreciable.  When  the  grid  becomes  dense  enough, 
the  discreteness  error  is  minimized,  as  shown  in  Fig.  A.4  (b).  For 
the  refined  grid,  the  convergent  simulation  results  of  FLUENT  and 
ANSYS  under  different  temperature  spans  are  shown  in  Fig.  A.5, 


where  numerical  results  of  the  performance  parameters  from  both 
models  display  almost  the  same  characteristics  in  this  ID  case. 

To  illustrate  how  the  surface  heat  loss  can  produce  multidimen¬ 
sional  effects,  the  3D  temperature  and  electric  potential  distribu¬ 
tions  in  the  thermocouple  by  the  FLUENT  model  are  depicted  in 
Figs.  A.6,  A.7  (a),  where  all  leg  surfaces  except  the  top  and  bottom 
junctions  are  assumed  to  be  exposed  to  heat  transfer  with  a  coef¬ 
ficient  of  500  W/m2  K,  representing  contributions  from  both  natu¬ 
ral  convection  and  radiation  heat  transfer.  The  bulk  mean 
temperature  in  the  module  cavity  is  set  to  be  the  arithmetic  aver¬ 
age  of  Th  and  Tc.  i.e.,  363  K.  With  the  same  model  parameters,  in 
Figs.  A.6,  A.7  (b)  3D  simulation  results  of  ANSYS  show  a  qualitative 
agreement  with  the  contours  from  the  FLUENT  model.  For  a  quan¬ 
titative  comparison,  the  performance  parameters  are  calculated  for 
both  models.  Due  to  the  prescribed  temperature  boundary  condi¬ 
tion  on  the  legs,  it  is  found  that  the  influence  of  surface  heat  loss 
on  P  and  /  is  negligible  in  this  case.  However,  with  the  surface  heat 


348 


363 


378 


393 


Fig.  A.6.  Simulation  of  temperature  (K)  distributions  for  a  convection  heat  transfer  coefficient  of  500  W/m2  K.  Th  and  Tc  are  423  K  and  303  K,  respectively,  (a)  FLUENT, 
(b)  ANSYS. 
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3  Contours  of  Potential  (V) 
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Fig.A.7.  Simulation  of  electric  potential  (addition  of  Ohm’s  and  Seebeck  potential,  V)  distributions  for  a  convection  heat  transfer  coefficient  of  500  W/m2  K.  Th  and  Tc  are  423  K 
and  303  K,  respectively.  The  cold  junction  of  the  p- type  (left)  leg  of  the  thermocouple  is  prescribed  as  ground  in  this  postprocessing,  and  a  scaled  load  (Rf  127)  is  assumed 
connected  between  the  cold  junctions  of  the  p-type  and  n-type  legs,  (a)  FLUENT,  (b)  ANSYS. 


loss,  Qh  is  significantly  increased  to  maintain  the  original  temper¬ 
ature  difference,  and  hence  r\  is  decreased.  The  changes  in  Qh  and 
ri  with  different  heat  transfer  coefficients  are  plotted  in  the 
close-up  of  Fig.  A.8,  where  both  numerical  models  display  almost 
the  same  characteristics  again  in  the  case  of  surface  heat  loss. 

For  the  device  example  modeled,  the  influence  of  the  non-ideal 
effects  of  the  thermal  and  electrical  interface  resistances  on  perfor¬ 
mance  parameters  turns  out  to  be  obvious  in  the  ID  study  [43]. 
These  interface  resistances  are  re-modeled  three-dimensionally 
in  conjunction  with  the  nonuniform  temperature  boundary  condi¬ 
tion  across  the  intermediate  components  of  the  soldering  bridges. 
The  detailed  profile  of  the  vertical  thermal  power  of  the  solid  mod¬ 
el  is  shown  in  Fig.  A.9.  The  thermal  and  electrical  interface  resis¬ 
tances,  and  the  Joule  heat  of  the  electrical  interface  resistance  are 
all  included  in  the  top  and  bottom  bridges,  although  the  last  has 
only  tiny  effects  in  this  case.  At  the  hot  and  cold  faces  of  the  ther¬ 
moelements  the  heat  conduction  has  a  sudden  change  in  the  verti¬ 
cal  direction,  reflecting  the  Peltier  heat  absorbed  and  evolved  at 
the  junctions  of  dissimilar  materials.  This  reversible  heat  is 


automatically  taken  into  account  by  the  thermoelectric  source 
term  described  in  Section  3. 

The  key  performance  parameters  calculated  by  the  model  with 
soldering  bridges  are  also  shown  in  Fig.  A.5  for  comparison. 
Although  there  are  still  some  differences  between  the  simulated 
and  measured  results  in  the  high  temperature  range  after  the  true 
3D  interface  is  modeled,  it  can  be  seen  in  Fig.  A.5  that  the  simula¬ 
tion  data  match  the  real  measurement  closely,  with  a  good  corre¬ 
lation.  The  acceptable  deviation  has  been  analyzed  and  can  be 
mainly  attributed  to  the  non-ideal  accuracy  of  the  heat  flow  deter¬ 
mination  in  the  device  test  rig  and  the  inherent  uncertainties  of  the 
commercial  thermoelectric  module  [43].  However,  the  negligible 
discrepancy  between  the  numerical  models  validates  the  proposed 
CFD  model  as  a  thermoelectric  simulator  equivalent  to  ANSYS,  and 
to  the  other  numerical  model  introduced  in  [43]  in  ID  cases.  Given 
that  the  experimental  determination  of  thermal  quantities  (such  as 
conductivity)  in  most  measurements  has  errors  of  a  similar  level, 
the  FLUENT  model  has  been  able  to  predict  the  efficiency  and  out¬ 
put  power  of  thermoelectric  generators  operating  in  fluid  systems 
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heat  transfer  coefficient  (W/m2K) 

Fig.  A.8.  Comparison  of  input  heat  flow  and  conversion  efficiency  between  FLUENT 
and  ANSYS  simulations  for  various  convection  heat  transfer  coefficients.  Th  and  Tc 
are  423  K  and  303  K,  respectively. 

with  sufficient  accuracy.  Further  refining  work  will  focus  on 
improving  the  accuracy  of  parameter  acquisition  in  the  test  rig 
for  the  model,  where  more  advanced  modules  will  be  measured 
under  hot  sources  of  higher  temperature. 

5.  Conclusions 

(i)  The  thermoelectric  processes  of  Seebeck,  Peltier,  and  Thom¬ 
son  effects  are  integrated  with  Joule  source  terms  through  a 
FVM  numerical  scheme  into  a  CFD  simulator.  The  proposed 
FLUENT  model  includes  the  temperature  dependence  of  all 
properties  of  the  p-  and  n-type  materials  composing  the 
thermocouple.  The  3D  modeling  results  provide  detailed 
profiles  of  temperature,  Seebeck  potential,  and  current  den¬ 
sity  as  well  as  the  values  of  power  and  efficiency.  Compari¬ 
sons  to  other  modeling  and  experimental  results  validate  the 
accuracy  of  the  numerical  model. 


(ii)  The  functionality  of  solving  scalar  fields  in  the  CFD  simulator 
has  been  extended  for  the  potential  field  in  solid  zones  of 
thermoelectric  generators.  The  power  source  model  is  com- 

^  prised  of  only  several  UDF  and  UDS,  and  hence,  it  is  scalable 

£  and  flexible  for  loading  in  FLUENT  to  interact  with  CFD  sub- 

models,  and  especially  useful  in  the  design  of  entire  power 
£=  systems  with  3D  temperature  and  heat  flux  profiles  [45]. 

^  The  numerous  existing  CFD  models  of  fluid  flow  and  com- 

2  bustion  in  FLUENT  can  be  immediately  connected  to  the 

B  thermoelectric  model  as  a  continuum  domain,  where  the 

cj 

a>  difficulties  encountered  in  the  implementation  of  the  multi- 

0  dimensional  boundary  condition  translation  of  the  fluid- 

structure  coupling  are  avoided.  To  handle  the  field-circuit 
interface,  a  general  flux  computation  is  defined  in  terms  of 
the  load  and  incorporated  into  the  UDS  boundary. 

(iii)  One  important  feature  of  the  present  model  is  the  incorpora¬ 
tion  of  the  reversible  contributions  from  inhomogeneous 
materials  and  the  Thomson  effect  into  the  source  term  of 
heat  conduction.  In  particular,  the  thermal  flux  in  the  com¬ 
putational  results  expresses  the  usual  conduction  heat  flow, 
and  thus,  it  has  a  direct  connection  with  the  temperature 
field  and  can  be  studied  separately  from  the  reversible  Pel- 
tier  and  Thomson  heat.  This  function  is  not  available  in 
ANSYS,  in  which  the  reversible  heat  and  the  irreversible  con¬ 
duction  are  inextricable  in  the  total  heat  flow.  Such  treat¬ 
ment  of  the  present  model  makes  relevant  simulation 
results  easier  to  analyze  and  to  truly  understand  as  com¬ 
pared  to  other  numerical  models  in  which  irreversible  and 
reversible  heat  are  defined  together  [37,40],  and  offers  a 
new  option  for  studying  the  effects  of  the  Thomson  heat. 

(iv)  Not  only  can  3D  calculations  in  the  thermoelectric  generator 
model  aid  in  optimizing  the  selection  of  the  3D  shape  and 
geometry  of  a  device,  the  effects  of  various  convection  and 
radiation  conditions  on  power  performance,  especially  in 
the  case  of  porous  medium,  can  also  be  easily  studied  by 
the  model  due  to  the  inherent  CFD  advantages  of  FLUENT 
in  modeling  and  implementing  such  source  terms.  Besides, 
the  3D  FLUENT  model  is  able  to  treat  both  isotropic  and 
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Fig.  A.9.  Simulation  of  vertical  heat  flux  (W/m2  in  Z  axis)  with  soldering  bridges.  Device  top  and  bottom  are  423  K  and  303  K,  respectively.  Peltier  heat’s  contribution  to  the 
interior  conduction  is  included  automatically  by  the  model. 
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orthotropic  physical  and  thermoelectric  properties  provided 
that  such  a  study  becomes  mandatory  to  the  overall  system 
behavior. 
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Appendix  A.  Supplementary  materials 

A  computational  case  and  the  UDF  library  associated  with  the 
proposed  model  can  be  found  in  the  on-line  version  as  the  supple¬ 
mentary  data  to  operate  in  Fluent  12  directly.  Supplementary  data 
associated  with  this  article  can  be  found,  in  the  online  version,  at 
doi:  1 0.1 01 6/j.ijheatmasstransfer.201 0.08.024. 
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